#ifndef BLOCKWISE_INTEGRATED_PRECALC_DIFFERENCES_HPP
#define BLOCKWISE_INTEGRATED_PRECALC_DIFFERENCES_HPP

namespace hms::swe
{


template<BlockLocation loc, int nRows,
	ReconstructionType recType, EdgeRange orien, typename LR_r>
auto phiDiff(
	ArrayXXs& phiDiffBuf,
	Index colBeg,
	const LR_r& cellLR, const LR_r& cellLLRR,
	const EdgeRow<recType, orien>& er
) -> std::array<ArrayNXsMap<nRows>, 3>
{
	using ER = EdgeRange;
	if constexpr ( orien == ER::horizontal ){
		assert( phiDiffBuf.size() >= er.nEdges * 12 );
		/* we calculate one block row of differences per row
		 * (except for the first row, when no differences are calculated yet),
		 * and cycle through them to reuse them:
		 * In the next iteration, er.i is incremented by 1, so the map pointers
		 * are moved by nRows buffer columns, therefore
		 * - the old diffM becomes the new diffL,
		 * - the old diffR becomes the new diffM and
		 * - the old diffL becomes invalid and the new diffL must be calculated.
		 * The differences must persist across block row iterations,
		 * therefore each variable must have its own memory space. */
		ArrayNXsMap<nRows>
		diffL {phiDiffBuf.col(colBeg + nRows*( (er.i + 0) % 3) ).data(),nRows,er.nEdges},
		diffM {phiDiffBuf.col(colBeg + nRows*( (er.i + 1) % 3) ).data(),nRows,er.nEdges},
		diffR {phiDiffBuf.col(colBeg + nRows*( (er.i + 2) % 3) ).data(),nRows,er.nEdges};
		assert(
			diffL.data() <= diffM.data() ?
			diffL.data() + diffL.size() <= diffM.data() : true );
		assert(
			diffM.data() <= diffR.data() ?
			diffM.data() + diffM.size() <= diffR.data() : true );
		assert(
			diffR.data() <= diffL.data() ?
			diffR.data() + diffR.size() <= diffL.data() : true );

		diffR = cellLLRR[1] - cellLR[1];
		/* during the first iteration, all differences must be calculated.
		 * For bottom blocks, this does not coincide with er.i == 0,
		 * as there only the last row of edges is reconstructed with a
		 * high order scheme. */
		if ( atBottom<loc>() || er.i == 0 ){
			diffL = cellLR[0] - cellLLRR[0];
			diffM = cellLR[1] - cellLR  [0];
		}
		#ifndef NDEBUG
		assert( diffL.isApprox(cellLR[0] - cellLLRR[0]) );
		assert( diffM.isApprox(cellLR[1] - cellLR  [0]) );
		#endif
		return {{
			std::move(diffL), std::move(diffM), std::move(diffR)
		}};
	} else if constexpr ( orien == ER::vertical ){
		assert(
			phiDiffBuf.data() + phiDiffBuf.size() >=
			phiDiffBuf.col(colBeg).data() + nRows*(2+er.nEdges)
		);
		/* for vertical edges, it's a lot simpler. The differences do not
		 * persist across block rows. We only have to subtract
		 * the cell values of the first cell until one before the end from
		 * the cell values of the second cell until the end.
		 * */
		ArrayNXsMap<nRows>
			diff  { phiDiffBuf.col(colBeg).data(), nRows, er.nEdges + 2 },
			diffL { phiDiffBuf.col(colBeg).data() + 0*nRows, nRows, er.nEdges },
			diffM { phiDiffBuf.col(colBeg).data() + 1*nRows, nRows, er.nEdges },
			diffR { phiDiffBuf.col(colBeg).data() + 2*nRows, nRows, er.nEdges };
		if constexpr ( atMiddleCols<loc>() ){
			assert(er.nEdges > 1); // otherwise the next line doesn't work
			diff.template leftCols<2>() =
				cellLR  [0].template leftCols<2>() -
				cellLLRR[0].template leftCols<2>();
			diff.rightCols(er.nEdges) = cellLLRR[1] - cellLR[1];
		} else {
			assert(er.nEdges == 1);
			diff.col(0) = cellLR  [0] - cellLLRR[0];
			diff.col(1) = cellLR  [1] - cellLR  [0];
			diff.col(2) = cellLLRR[1] - cellLR  [1];
		}
		return {{
			std::move(diffL), std::move(diffM), std::move(diffR)
		}};
	}
}

} // namespace hms::swe

#endif
